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Abstract 

Based  on  the  notion  Minimizing  Lipschitz  Extensions  and 
its  connection  with  the  infinity  Laplacian,  a  computational 
framework  for  surface  warping  and  in  particular  brain 
warping  (the  nonlinear  registration  of  brain  imaging  data) 
is  presented  in  this  paper.  The  basic  concept  is  to  com¬ 
pute  a  map  between  surfaces  that  minimizes  a  distortion 
measure  based  on  geodesic  distances  while  respecting  the 
boundary  conditions  provided.  In  particular,  the  global 
Lipschitz  constant  of  the  map  is  minimized.  This  frame¬ 
work  allows  generic  boundary  conditions  to  be  applied 
and  allows  direct  surface-to-surface  warping.  It  avoids  the 
need  for  intermediate  maps  that  flatten  the  surface  onto 
the  plane  or  sphere,  as  is  commonly  done  in  the  litera¬ 
ture  on  surface-based  non-rigid  brain  image  registration. 
The  presentation  of  the  framework  is  complemented  with 
examples  on  synthetic  geometric  phantoms  and  cortical 
surfaces  extracted  from  human  brain  MRI  scans. 

1  Introduction 

Brain  warping,  a  form  of  brain  image  registration  and 
geometric  pattern  matching,  is  one  of  the  most  funda¬ 
mental  and  thereby  most  studied  problems  in  computa¬ 
tional  brain  imaging  [37].  Brain  images  are  commonly 
warped,  using  3 -dimensional  deformation  fields,  onto  a 
common  neuroanatomic  template  prior  to  cross-subject 
comparison  and  integration  of  functional  and  anatomi¬ 
cal  data.  Images  of  the  same  subject  may  be  warped 
into  correspondence  over  time,  to  help  analyze  shape 
changes  during  development  or  degenerative  diseases.  Al¬ 
most  all  the  active  research  groups  in  this  area  have  de¬ 
veloped  and/or  have  their  favorite  brain  warping  tech¬ 


nique.1  A  few  representative  works  can  be  found  at 
[5,  7,  8,  13,  14,  16,  30,  34,  36,  37,  40,  41],  this  list  be¬ 
ing  far  from  complete.  In  spite  of  this,  the  problem  is 
still  open  and  widely  studied,  since  there  is  not  a  “ground 
truth”  method  to  obtain  a  map  between  brains.  The  cri¬ 
teria  for  matching  different  features  (e.g.,  geometry  or  in¬ 
tensity)  may  also  depend  on  the  applications,  which  range 
from  recovering  intraoperative  brain  change  to  mapping 
brain  growth,  or  reducing  cross-subject  anatomical  differ¬ 
ences  in  group  functional  MRI  studies. 

The  way  the  brain  warping  problem  is  addressed  is 
critical  for  studies  of  brain  diseases  that  are  based  on 
population  comparisons.  Examples  of  this  application 
can  be  found  at  [14,  34],  although  these  are  a  very  non- 
exhaustive  account  of  the  rich  literature  on  the  subject. 
The  interested  reader  may  also  check  [35]  for  numer¬ 
ous  applications  of  brain  warping  and  population  stud¬ 
ies.  As  detailed  in  [37],  brain  warping  approaches  can 
be  divided  into  two  classes,  those  based  on  volume-to- 
volume  matching  and  those  based  on  surface-to-surface 
matching.  Our  work  belongs  to  the  latter  of  the  cate¬ 
gories.  Surface  matching  has  recently  received  increas¬ 
ing  attention  as  most  functional  brain  imaging  studies  fo¬ 
cus  on  the  cortex,  which  varies  widely  in  geometry  across 
subjects.  The  power  of  these  studies  depends  on  the  de¬ 
gree  to  which  the  functional  anatomy  of  the  cortex  can  be 
aligned  across  subjects,  so  improved  cortical  surface  reg¬ 
istration  has  become  a  major  goal.  In  contrast  with  flow 
based  works  such  as  those  in  [5,  30,  34],  our  motivation 
is  as  in  [1,  16,  17,  18,  21,  38,  39,  42].  That  is,  we  aim  to 
compute  a  map  that  preserves  certain  pre-defined  geomet¬ 
ric  characteristics  of  the  surfaces.  While  the  literature  has 


lrThis  includes  groups  at  JHU,  UCLA,  U.  Penn.,  INRIA,  MGH,  GAT- 
ECH,  Harvard-B  W,  and  the  University  of  Florida,  to  name  just  a  few. 


mainly  attempted  to  preserve  angles  and  areas,  we  work 
with  geodesic  distances  (see  also  [33]).  Our  work  is  in¬ 
spired  by  the  literature  on  Lipschitz  minimizing  maps  and 
in  its  connection  to  the  infinite  Laplacian.  The  motivation 
for  using  these  frameworks  will  be  presented  after  some 
brief  mathematical  introduction  below. 

In  this  paper  we  therefore  introduce  the  use  of  Lip¬ 
schitz  minimizing  maps  into  the  area  of  computational 
brain  imaging,  presenting  a  theoretical  and  computational 
framework  complemented  by  examples  with  artificial  and 
real  data.  An  additional  critical  contribution  of  the  work 
here  presented  is  that  intermediate  distorting  maps  to  the 
plane  or  sphere  are  avoided  -  these  intermediate  mappings 
are  common  practice  in  the  brain  warping  literature. 

2  Formal  statement  of  the  problem 

Let  Bi  and  S2  be  two  cortical  surfaces  (2D  surfaces  in  the 
three  dimensional  Euclidean  space)  which  we  consider 
smooth  and  endowed  with  the  metric  inherited  from  1R3 
so  that  dg1  and  dg2  are  the  geodesic  distances  measured 
on  B\  and  B2 ,  respectively.  Let  Ti  C  B\  and  T2  C  B2  be 
subsets  which  represent  features  for  which  a  correspon¬ 
dence  is  already  known.  In  general,  the  sets  Ti  are  the 
union  of  smooth  curves  traced  on  the  surfaces,  e.g.,  sulcal 
beds  lying  between  gyri,  and/or  a  union  of  isolated  points. 
A  set  of  anatomical  landmarks  that  occur  consistently  in 
all  subjects  can  be  reliably  identified  using  standardized 
anatomical  protocols  or  automated  sulcal  labeling  tech¬ 
niques  (see  for  example  Brain  VISA  by  Mangin  and  Riv¬ 
iere  and  SEAL  by  Le  Goualher). 

Functional  anatomy  also  varies  with  respect  to  sulcal 
landmarks,  but  sulci  typically  lie  at  the  interfaces  of  func¬ 
tionally  different  cortical  regions  so  aligning  them  im¬ 
proves  the  registration  of  functionally  homologous  areas. 
As  commonly  done  in  brain  warping  [37],  we  assume  that 
a  correspondence  between  T 1  and  T2  is  pre-specified  to 
the  map  (boundary  conditions  of  the  map).  In  this  cor¬ 
respondence,  internal  point  correspondences  may  be  al¬ 
lowed  to  relax  along  landmark  curves  in  the  final  map¬ 
pings,  e.g.,  [26]. 

To  fix  ideas  let’s  assume  that  T 1  =  U k=ixi  and  T2  = 
UfLi  Vi’  and  that  the  correspondence  is  given  by  X{  yi 
for  1  <  i  <  N. 

We  want  to  find  a  (at  least  continuous)  map  f  :  B\  -A 


B2  such  that  (p(xi)  —  yi  for  1  <  i  <  N  and  such  that  f 
produces  minimal  distortion  according  to  some  functional 
J.  One  possible  way  of  interpreting  this  problem  is  that 
we  are  trying  to  extrapolate  or  extend  the  correspondence 
from  T 1  to  the  whole  of  B\  in  such  a  way  that  we  achieve 
small  distortion. 

A  possible  way  to  measure  the  distortion  produced  by 
a  map  f  is  by  computing  the  functionals  (1  <  p  <  00) 


where  Djg1  denotes  differentiation  intrinsic  to  the  surface 
B\  and  p  is  the  area  measure  on  B\ .  One  immediate  idea 
is  then  to  consider,  for  a  fixed  p  G  (1,  00),  the  following 
variational  problem: 

Problem  1  (minimize  Jp)  Find  G  S  such  that 
J p((j))  =  inf^es  where  S  is  a  certain  smoothness 

class  of  maps  f  from  B\  to  B2  such  that  they  respect  the 
given  boundary  conditions  f(xi)  =  yi  for  all  X{  G  Ti. 

The  case  p  —  2  corresponds  to  the  Dirichlet  functional 
and  has  connections  with  the  theory  of  (standard)  Har¬ 
monic  Maps.  In  more  generality,  it  is  customary  to  call 
the  solutions  to  Problem  1  p-Harmonic  Maps,  see  for  ex¬ 
ample  [10,  20,  11].  It  is  easy  to  show,  under  mild  regular¬ 
ity  assumptions,  that  for  a  fixed  Jp((p)  is  nondecreasing 
as  a  function  of  p,  and  that  [15] 

J00O)  :=  lim  Jp(4>)  =  essupxeB  \\DBl4>(x)\\2,  (2) 

pjoo 

which  is  the  Lipschitz  constant  of  <fi. 

In  this  paper  we  propose  to  use  the  functional  Jqo  as  a 
measure  of  distortion  for  maps  between  cortical  surfaces 
and  to  solve  the  associated  variational  problem  in  order 
to  find  a  candidate  mapping  between  the  cortical  surfaces 
(constrained  by  the  provided  boundary  conditions). 

Let  C  denote  the  space  of  all  Lipschitz  continuous  maps 
ft  :  B\  — y  B2  such  that  ^{xf)  =  yi  for  1  <  i  <  N.  We 
then  propose  to  solve  the  following  problem: 

Problem  2  (minimize  Jqo)  Find  G  C  such  that 

Joo(0)  =  ifrfpEvC  Joo 

We  now  argue  in  favor  of  this  functional. 
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2.1  Why  use  Joo  ? 

Our  first  argument  is  that  measures  distortion  in  a 
more  global  way  than  any  of  the  Jp  for  p  E  (l,oo), 
since  instead  of  computing  an  averaged  integral  quantity, 
we  are  looking  at  the  supremum  of  the  local  distortions, 
\\dBi4>{x)\\2.  Note  also  that  as  stated  above,  upper- 
bounds  Jp  under  mild  regularity  assumptions. 

Another  element  to  consider  is  that  this  problem  is  well 
posed  for  the  kind  of  general  boundary  data  we  want  to 
respect,  provided  both  at  curves  and  isolated  points  on  the 
surfaces.  At  least  for  the  case  p  <  2,  this  is  not  true  in 
general,  see  [3]. 

We  are  then  looking  for  a  Lipschitz  extension  of  the 
map  given  at  Ti  whose  Lipschitz  constant  is  as  small  as 
possible.  Let 


i(ri,r2) 


Xi,Xj£zY\  Cl]Si  \%ii  % j  ) 


that  is,  the  Lipschitz  constant  of  the  boundary  data.  In 
general,  we  will  have  inf^€£  >  L(ri,r2).  This 
is  related  to  Kirszbraun’s  Theorem,  which  in  one  of  its 
many  guises  states  that  a  Lipschitz  map  f  :  S  MD , 

S  C  Md ,  has  an  extension  /  :  Md  1RD  with  the 
same  Lipschitz  constant  as  /,  see  [12].  In  the  same  vein, 
one  has  Whitney  and  McShane  extensions  which  apply 
to  the  case  when  the  domain  is  any  metric  space  X  and 
the  target  is  1R.  These  extensions  provide  functions  that 
agree  with  /  where  boundary  conditions  are  given  and 
preserve  the  Lipschitz  constant  throughout  X,  see  for  ex¬ 
ample  [2,  22].  The  more  general  problem  of  extending 
/  :  S  4  F  (S  C  X,  X  and  Y  any  metric  spaces)  to  all 
X  with  the  same  Lipschitz  constant  is  not  so  well  under¬ 
stood  and  only  partial  results  are  known,  see  for  example 
[23,24,25]. 

The  idea  then  is  to  keep  the  distortion  at  the  same  order 
as  that  of  the  provided  boundary  conditions.  In  general 
there  might  be  many  solutions  for  the  Problem  (2).  One 
particular  class  of  minimizers  which  has  recently  received 
a  lot  of  attention  is  that  of  absolute  minimizers,  or  abso¬ 
lutely  minimizing  Lipschitz  extensions  (AMLE).  Roughly 
speaking,  the  idea  here  is  to  single  out  those  solutions 
of  Problem  (2)  that  also  possess  minimal  local  Lipschitz 
constant,  again,  see  [2]  for  a  general  exposition,  and  [22] 
for  a  treatment  of  the  case  when  the  domain  is  any  reason¬ 
able  metric  space  and  the  target  is  the  real  line. 


3  Proposed  computational  ap¬ 
proach 

If  we  take  for  example  the  case  of  p-Harmonic  maps,  one 
way  of  dealing  with  the  computation  of  the  optimal  map 
4>p  is  by  implementing  the  geometric  p-heat  flow  associ¬ 
ated  with  the  Euler-Lagrange  equation  of  the  functional 
Jp,  starting  from  a  certain  initial  condition.  As  was  ex¬ 
plained  in  [29],  using  an  implicit  representation  for  both 
B\  and  S2,  we  could  obtain  the  partial  differential  equa¬ 
tion  PDEp  we  need  to  solve  in  order  to  find  op.  By  tak¬ 
ing  the  formal  limit  as  p  |  oo  we  would  find  PDEqo, 
the  PDE  that  characterizes  the  solution  ^oo  of  the  (varia¬ 
tional)  Problem  (2). 2  All  of  this  might  work  if  we  had  a 
notion  of  solution  for  the  resulting  PDEs.  Whereas  this 
is  feasible  in  the  case  of  PDEp  for  1  <  p  <  oo,  to  the 
best  of  our  knowledge,  there  is  no  such  notion  of  a  so¬ 
lution  for  PDEqo.  One  could  of  course  still  persist  and 
try  to  solve  these  equations  without  the  necessary  theo¬ 
retical  foundations  and  call  these  plausible  solutions  oo- 
Harmonic  Maps.  Nonetheless,  this  is  certainly  an  inter¬ 
esting  line  of  research. 

A  different  direction  is  considered  in  this  work.  As  a 
guiding  example,  we  first  concentrate  on  the  case  where 
B\  is  any  closed  smooth  manifold  and  B2  is  replaced  by 
M,  as  considered  in  [4]  (for  scalar  data  interpolation  on 
surfaces),  and  in  [32].  In  [4],  the  authors  propose  to 
follow  a  similar  path  to  the  one  we  have  just  described, 
and  they  do  not  obtain  a  convergent  numerical  discretiza¬ 
tion  for  the  resulting  PDE.  Meanwhile,  in  [32],  the  author 
proposes  a  convergent  discretization  of  the  PDE,  basing 
his  construction  on  the  original  variational  problem.  We 
choose  to  follow  this  idea  as  our  guiding  principle. 

We  now  explain  this  alternative  approach.  The  ba¬ 
sic  idea  is  simple,  instead  of  first  obtaining  the  Euler- 
Lagrange  equations  for  and  then  discretizing  them, 
we  will  first  discretize  and  then  proceed  to  solve 
the  resulting  discrete  problem.  Consider  that  the  do¬ 
main  B\  is  given  discretely  as  a  set  of  (different)  points 
Bi  =  {xi, . . .  jXm}  together  with  a  neighborhood  rela¬ 
tion  (i.e.,  a  graph).  To  fix  ideas  let’s  assume  the  neighbor¬ 
hood  relation  is  a  /^-nearest  neighbors  one.  Denote,  for 
each  1  <  i  <  m,  by  Ni  =  {xj1 , . . . ,  Xjk  }  G  li  the  set 

2The  case  when  the  domain  is  a  subset  of  lRd  and  the  target  is  the 
real  line  leads  to  the  so  called  infinity  Laplacian,  see  [2,  19,  9]. 
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of  k  neighbors  of  the  point  X{ .  We  consider  the  discrete 
local  Lipschitz  constant  of  the  map  0  at  Xi : 


U{<p) 


max 

x3eNi 


dBMKXzJMXj)) 

dBl  (xi,xj) 


(3) 


Upon  noting  that  Li  (0)  serves  as  a  discrete  approxima¬ 
tion  to  1 1  Db1  (j)(xi )  1 1 2 ,  we  see  that  a  possible  discretization 
of  the  functional  Joo(^)  is  given  by  the  discrete  global 
Lipschitz  constant  of  0  given  by  maxi<^<m  Li  ((/>).  The 
author  of  [32]  proposed,  in  the  case  when  S2  is  replaced 
by  1R ,  solving  the  discrete  version  of  Problem  (2)  by  fol¬ 
lowing  the  following  iterative  procedure  (here  described 
for  S2  a  surface  as  in  our  problem): 


point.  For  computational  efficiency,  we  work  at  all  times 
with  two  different  scales  in  the  discrete  domain  Bi .  We 
choose  a  subset  F\  of  Bi  such  that  <C  m  but  still 
Fi  is  an  efficient  (well  separated)  covering  of  Bi  with 
small  covering  radius.  We  do  this  by  using  the  well  known 
(geodesic)  Farthest  Point  Sampling  (FPS)  procedure,  see 
[28,  31],  which  can  be  efficiently  constructed  based  on 
optimal  computational  techniques.  Roughly  speaking,  we 
apply  the  iterative  procedure  on  this  subset  of  points  only 
and  then  extend  the  map  to  the  rest  of  the  points  in  the  do¬ 
main  Bi .  We  now  show  how  to  obtain  a  reasonable  initial 
condition  (j)0  and  then  discuss  additional  details  regarding 
the  implementation  of  the  iterative  procedure  described  in 
the  previous  section. 


•  Let  4> 0  be  an  initial  guess  of  the  map. 


•  For  each  n  >  1,  if  Xi  £  Ti,  let 


4>n(%i)  =  argminmax 
yeJ32 


dB2(y,<t>n-i{xj)) 

dBl  (• Xi,Xj ) 


•  c t>n{xi )  =  Vi  for  all  n  >  0  for  Xi  eT  1. 


Building  the  initial  condition:  We  compute,  for  all  xr  G 

^i\ri>  <f>o(xr)  =  &vgmmy€B2  max^er,  '  For 

this  step  we  use  the  classical  Dijkstra’s  algorithm  for  ap- 
(4)  proximating  the  distances  d^1  and  d&2  since  they  might 
be  evaluated  at  faraway  points.  This  is  of  course  run  on 
the  graphs  obtained  from  connecting  each  point  to  its  k- 
nearest  neighbors. 


With  computational  efficiency  related  modifications 
described  below,  this  is  the  approach  we  follow  in  gen¬ 
eral.  The  intuition  behind  this  iterative  procedure  is  that, 
at  each  point  of  the  domain,  we  are  changing  the  value  of 
the  map  in  order  to  minimize  the  local  Lipschitz  constant, 
that  is,  the  local  distortion  produced  by  the  map.  This  is  in 
agreement  with  the  notion  of  AMLEs  briefly  explained  in 
§2.1.  We  should  remark  that  since  we  are  using  intrinsic 
distances  for  the  matching,  we  can  let  Li{(j))  play  the  role 
of  (the  norm  of)  the  displacement  field  for  analyzing  the 
deformation,3  see  §5  ahead. 

4  Implementation  details 


The  iterative  procedure:  After  is  computed  for  all 
points  in  the  set  F\ ,  we  run  the  iterative  procedure  from 
§2  on  this  set  of  points.  The  main  modification  here  is  that 
whereas  we  still  use  Dijkstra’s  algorithm  for  approximat¬ 
ing  d&2  in  the  target  surface,  since  in  the  domain  we  must 
compute  d&1  only  for  neighboring  points  (T\  was  chosen 
to  be  dense  enough),  for  computationally  efficiency  we 
can  approximate  dg1(xi,Xj)  ~  || X{  —  x3  \ |  for  x3  E  N%. 
We  should  also  point  out  that  for  points  in  F\ ,  the  neigh¬ 
borhood  relation  is  defined  to  be  that  of  /^-nearest  neigh¬ 
bors  with  respect  to  the  metric  on  Bi  defined  by  the  adja¬ 
cency  matrix  of  Bi .  Let  :  F\  B2  denote  the  map 
obtained  as  the  output  of  this  stage. 


In  addition  to  discretizing  the  domain  Bi ,  we  also  use  a 
discretization  B2  =  {y± , . . . ,  ymi }  of  the  target  space  B2 
for  our  implementation.  We  endow  B2  with  a  neighbor¬ 
hood  relation  given  by  the  /^-nearest  neighbors  of  each 

3  One  can  imagine  a  situation  in  which  two  isometric  surfaces  are 
matched  by  our  algorithm  such  that  Li(<f>)  =  1  for  all  i,  but  the  dis¬ 
placement  held  || Xi  —  (ft(xi)\  \  is  large  since  there  may  be  no  rigid  mo¬ 
tion  that  aligns  the  two  surfaces.  One  simple  example  is  a  hat  sheet  of 
paper  and  the  same  sheet  slightly  bent. 


Extension  to  the  whole  domain:  After  we  have  iterated 
over  points  in  F\  until  convergence,  we  extend  the  map 
<; 1 L  to  all  points  X{  in  Bi  \{Fi  U  Ti }.  This  is  done  by  com- 
puting  4>*{Xi)  =  argmin^^maxseFiur!  £Bl(XitX)  • 
For  this  step,  and  since  we  have  already  obtained  the  map 
for  a  relatively  dense  subset,  we  approximate  both  d&1 
and  ds2  by  the  Euclidean  distance.  Once  again,  the  moti¬ 
vation  for  this  is  just  computational  efficiency. 
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5  Examples 

In  this  section  we  present  some  computational  examples 
of  the  ideas  presented  in  previous  sections.  First,  in  Fig¬ 
ure  1  the  domain  B\  is  a  cube  ( m  =  10086)  and  the  target 
B 2  is  a  sphere  (mr  =  17982).  For  the  purposes  of  visual¬ 
izing  the  map,  we  assigned  the  clown  texture  (which  can 
be  thought  of  as  a  function  I  :  B2  M)  to  the  sphere, 
which  can  be  seen  on  the  bottom-right  corner  of  the  fig¬ 
ure.  The  sphere  and  the  cube  were  concentric  and  of  ap¬ 
proximately  the  same  size.  We  selected  F\  on  the  cube 
consisting  of  1000  well  separated  points  using  the  FPS 
procedure  alluded  to  in  §4.  Also,  we  set  k  =  6  (num¬ 
ber  of  neighbors).  We  then  chose  T i  to  be  the  first  100 
points  of  the  set  and  then  projected  them  onto  the  sphere, 
obtaining  in  this  way,  the  corresponding  set  T2  to  use  as 
boundary  conditions.  We  then  followed  the  computational 
procedure  detailed  before.  The  top-left  figure  shows  the 
composition  I  o  :  cube  — ^  M  as  a  texture  on  the  cube. 
Finally,  the  top-right  and  the  bottom-left  images  show  the 
histogram  of  Li((j) *)  and  its  spatial  distribution  in  the  do¬ 
main  (we  paint  the  cube  at  each  point  Xi  with  the  color 
corresponding  to  Li(q !>*)),  respectively.  Ideally,  we  would 
like  to  obtain  a  (5-type  histogram,  meaning  that  the  dis¬ 
tances  have  been  constantly  scaled.  Of  course,  this  is  not 
possible  (unless  one  of  the  surfaces  is  isometric  to  a  scaled 
version  of  the  other),  and  we  attempt  to  obtain  histograms 
as  concentrated  as  possible.  This  is  quite  nicely  obtained 
for  this  and  the  additional  examples  in  this  paper. 

Figure  2  shows  the  construction  of  a  map  from  the 
unit  sphere  S 2  into  a  cortical  hemisphere  B  (£>).  The 
boundary  conditions  consisted  of  6  pairs  of  points.  We 
first  took  the  following  6  points  on  the  sphere  T  i  = 
{(=bl,  0, 0),  (0,  ±1,  0),  (0, 0,  ±1)}.  We  then  constructed 
the  intrinsic  distance  matrix  [ds^  (pi,Pj)\  for  all  pl ,  pj  E 
T i.  Finally,  we  chose  6  points  {gi, . . .  ,qQ}  =  T2  in  B 
such  that  max^j  was  as  c^ose  as  Possible  to 

^diam(B).  We  painted  B  with  a  texture  In  depending 
on  its  mean  curvature  so  as  to  more  easily  visualize  the 
sulci/crests:  If  H (x)  stands  for  mean  curvature  of  B  at  x, 
then  Ih(x)  =  (H(x)  —  min^  H(x))2 .  See  the  caption 
for  more  details. 

The  example  in  Figure  3  is  about  computing  a  map  4> 
from  a  subject’s  left  hemisphere  Bi  to  another  subject’s 
left  hemisphere  B2 .  The  boundary  conditions  were  con¬ 


structed  in  a  way  similar  to  the  one  used  for  the  previ¬ 
ous  example,  but  in  this  case,  300  points  were  chosen. 
Note  that,  if  available,  hand  traced  curves  could  be  used 
as  commonly  done  in  the  literature  (in  other  words,  more 
anatomical/functional  oriented  boundary  conditions).  In 
the  first  two  rows  we  show  4  different  views  of  each  corti¬ 
cal  surface,  and  in  the  third  row  we  show  Bi  colored  with 
the  values  of  Ll(&)  which  we  interpret  as  a  measure  of 
the  local  deformation  of  the  map  needed  to  match  Bi  to 
B2 .  See  the  caption  for  more  details. 


Figure  1 :  Artificial  example  of  the  proposed  warping  algo¬ 
rithm.  From  top  to  bottom  and  left  to  right:  The  domain  surface, 
with  a  picture  painted  on  it  to  help  in  visualizing  the  computed 
map;  histogram  of  the  Lipschitz  constant  ( note  how  it  is  con¬ 
centrated  around  a  single  value);  color  coded  distribution  of  the 
Lipschitz  constant  for  the  computed  map;  and  mapped  texture 
following  the  computed  map. 


6  Concluding  remarks 

In  this  paper  we  have  introduced  the  notions  of  minimiz¬ 
ing  Lipschitz  extensions  into  the  area  of  surface  and  brain 
warping.  These  maps  provide  a  more  global  constraint 
than  ordinary  p-harmonic  ones,  and  allow  for  more  gen¬ 
eral  boundary  conditions.  The  proposed  computational 
framework  leads  to  an  efficient  surface-to-surface  warp¬ 
ing  algorithm  that  avoids  distorting  intermediate  steps  that 
are  common  in  the  brain  warping  literature.  We  are  cur¬ 
rently  investigating  the  use  of  this  new  warping  technique 
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for  creating  population  averages  and  applying  it  to  dis¬ 
ease  and  growth  studies.  In  earlier  work,  the  Jacobian  of 
a  deformation  mapping  over  time  has  been  used  to  map 
the  profile  of  brain  tissue  growth  and  loss  in  a  subject 
scanned  serially  (tensor-based  morphometry  [6,  34]).  The 
discrete  local  Lipschitz  constants  of  our  computed  map¬ 
pings  also  provide  a  useful  index  of  deformation  that  can 
be  analyzed  statistically  across  subjects.  The  framework 
here  introduced  can  also  be  applied  in  3D  for  volumetric 
warping  and  with  weighted  geodesic  distances  instead  of 
natural  ones  to  include  additional  geometric  characteris¬ 
tics  in  the  matching.  As  recently  shown  in  [27],  the  use  of 
pairwise  distances  is  of  importance  of  other  matching  and 
computer  vision  tasks.  Results  in  these  directions  will  be 
reported  elsewhere. 
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Figure  2:  An  example  of  mapping  between  the  cortex  and  a 
sphere.  The  order  is  the  same  as  in  the  previous  figure,  but 
now  the  domain  and  target  surface  are  colored  with  a  curvature- 
based  color  code.  Note  once  again  the  concentration  of  the  Lip- 
schitz  constant  for  the  computed  map.  On  the  top,  the  texture 
map  corresponding  to  Ih(x)  as  described  in  the  text  is  used. 
On  the  bottom,  the  texture  is  ma x(I h  (x) ,  5)  for  a  user  selected 
value  of  the  threshold  5,  which  highlights  the  gyral  crests. 
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Figure  3:  Warping  between  the  cortical  surfaces  of  two  brains.  In  the  first  row  we  show  4  views  of  Bi :  posterior,  medial,  lateral 
and  directly  viewing  the  occipital  cortex.  The  corresponding  4  views  of  B2  are  shown  in  the  second  row.  In  the  third  row,  we 
show  Bi  with  texture  I(xi)  =  Lj($)  which  can  interpreted  as  a  measure  of  local  deformation  needed  to  match  Xi  G  Bi  to 
<&(xi)  G  B2.  Relatively  little  deformation  (blue  colors)  is  required  to  match  features  across  subjects  on  the  flat  interhemispheric 
surface  ( second  image  in  the  second  row).  This  is  consistent  with  the  lower  variability  of  the  gyral  pattern  in  the  cingulate  and 
medial  frontal  cortices.  By  contrast,  there  is  significant  expansion  required  to  match  the  posterior  occipital  cortices  of  these  two 
subjects,  especially  in  the  occipital  poles  which  are  the  target  of  many  functional  imaging  studies  of  vision.  The  final  panel  in  the 
figure  shows  the  corresponding  histogram  for  Li{(jj),  the  local  Lipschitz  constants  of  the  map. 
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